clear all;

%% declare variables
filename = 'e3registro.txt';
report = 'report.txt';
clases = 7;

for i=1:clases
	cant(i) = 0;
	O(i) = 0;
	M(i) = 0;
	X(i) = 0;
end%for

%%abro archivo para el report
report = fopen(report, 'w');
fprintf(report, 'Resultados para E3\n\n');

%% abro, leo el archivo y lo cierro
file = fopen(filename);
fscanf(file, '%c', 66);	%%consumo el header
i = 1;
time_max = 0;
time_min = 1000;
%while (!feof(file)) % octave ...
while (~ feof(file))
	times(i) = fscanf(file, '%f', 1);
	if (times(i) > time_max)
		time_max = times(i);
    end%if
	if (times(i) < time_min)
		time_min = times(i);
    end%if
	%i++;
    i = i + 1;
end%while
fclose(file);
n = i-1;


%% computo las clases
dif = (time_max - time_min);
fprintf(report, 'N\t\t%d\n', n);
fprintf(report, 'Clases\t\t%d\n', clases);
fprintf(report, 'Maximo\t\t%g\nMinimo\t\t%g\n', time_max, time_min);
for i=1:n
	for j=1:clases
		if (times(i) > ((j-1)/clases)*dif+time_min && times(i) < (j/clases)*dif+time_min)
			O(j) = O(j) + 1;
        end%if
		M(j) = ((j-1)+0.5)*((time_max-time_min)/clases) + time_min;
    end%for
end%for


%% computo los estadisticos de maxima verosimilitud
mu = 0;
sigma = 0;
for j=1:clases
	X(j) = M(j) * O(j);
	mu = mu + X(j);
end%for
mu = mu / n;
for j=1:clases
	sigma = sigma + (X(j) - mu)^2;
end%for
sigma = sqrt(sigma/n);

sigma = 0;
for i=1:n
	sigma = sigma + ((times(i) - mu)^2);
end%for
sigma = sqrt(sigma/n);
fprintf(report, 'mu\t\t%g\nsigma\t%g\n', mu, sigma);


%% calculo la cantidad esperada para la normal
fprintf('j\tO(j)\tM(j)\t\tX(j)\t\ta\t\tb\t\tE(j)\n');
aux=zeros(1,clases);
for j=1:clases
    X(j) = O(j)*M(j);
    a = ((j-1)/clases)*dif+time_min;
    b = (j/clases)*dif+time_min;
    aux(j)=b;
    E(j) = (normcdf(b,mu,sigma) - normcdf(a,mu,sigma)) * n;
    fprintf('%d\t%d\t%f\t%f\t%f\t%f\t%f\n',j,O(j),M(j),X(j),a,b,E(j));
end%for



%% TEST CHI CUADRADO
%% H0: PROVIENE DE UNA EXPONENCIAL
alpha = 0.05   % probability that Chi^2 exceeds
               % critical value c = P(Chi^2 >= c).
SumChi = 0;
fprintf(report, '\nTEST CHI CUADRADO\n');
fprintf(report, 'H0: El tiempo entre arribos de E3 proviene de una distribucion normal N ~ (mu, sigma)\n');
fprintf(report, 'Clases\tO(i)\tE(i)\t\tO(i)-E(i)\t(O(i)-E(i))^2\t\t(O(i)-E(i))^2 / E(i)\n');
for j=1:clases
	DifCuad(j) = ((O(j)-E(j))^2);
	Chi(j) = DifCuad(j) / E(j);
	SumChi = SumChi + Chi(j);
    fprintf(report, '%d\t%d\t%g\t\t%g\t\t%g\t\t\t%g\n',j,O(j),E(j),O(j)-E(j),DifCuad(j),Chi(j));
end%for
fprintf(report, '\nEl estadistico Chi Cuadrado es %g\n',SumChi);
fprintf(report, 'Buscar en tabla el Chi Critico para %d grados de libertad y el nivel de confianza deseado\n', clases-1);

df = clases - 2 % degree of freedom ...
c = getCriticalValue_Chi2(alpha, df);
checkNullHypothesis(SumChi, c);

% check the result with the matlab-function "chi2gof" (to compare if we get the same result) ...
% which solution is correct?
O
E
[h, p, stats] = chi2gof(aux, 'ctrs', aux, 'frequency', O, ...
                             'expected', E, 'nparams', 1)
c = getCriticalValue_Chi2(alpha, stats.df);
checkNullHypothesis(stats.chi2stat, c);

% --> this calculates the correct values!
[h2, p2, stats2] = chi2gof(times, 'cdf', {@normcdf, mean(times), std(times, 1)}, 'nbins', clases)
% [h2, p2, stats2] = chi2gof(times, 'cdf', @(x) normcdf(x, mean(times), std(times, 1)), ...
%                                   'nbins', clases, 'nparams', 2) % the same (optional) ...
c = getCriticalValue_Chi2(alpha, stats2.df);
checkNullHypothesis(stats2.chi2stat, c);


						  
%% TEST Q-Q
y = sort(times);
for i=1:n
	gamma(i) = (i-0.5)/n;
	gammaInv(i) = norminv(gamma(i), mu, sigma);
end%for

%% create the Q-Q plot:
iStdSize = 12;

figure(1);
clf;
hold on;

plot(gammaInv, y, 'ko', 'LineWidth', 1, 'MarkerEdgeColor', 'k', ...
                  'MarkerFaceColor', 'w', 'MarkerSize', 5);
n = 1.7:0.01:2.4;
plot(n, n, '-k', 'LineWidth', 1);

%legend('Plot Q-Q');
%xlabel('Quantil Teorico');
%ylabel('Quantil Observado');
%print -deps -mono'PlotQQ_E3.eps';

xlbl = xlabel('\it Quantil Teorico');
ylbl = ylabel('\it Quantil Observado');
set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
axis equal;
%axis([1.7 1.7 2.4 2.4]);
grid on;

hold off;

%% histograma
figure(2);
clf;
hold on;

%histfit(times, clases);
hist(times, clases);
%xRange = floor(min(times)) : 0.01 : ceil(max(times));
% xRange = linspace(1.7, 2.3, 100);
% pdf_data = clases*normpdf(xRange, mean(times), std(times, 1));
% plot(xRange, pdf_data, 'k--');

%pdf_data2 = clases*normpdf(aux, mu, sigma);
%plot(aux, pdf_data2, 'k.');

h = findobj(gca,'Type','patch');
set(h,'FaceColor', 'w','EdgeColor','k')
l = legend( '\it Histograma de tiempos entre llegadas al registro');%, '\it PDF');
set( l, 'Interpreter', 'tex', 'Location', 'NorthWest', ...
        'FontName', 'Times', 'FontSize', 9 );
% correct the box size of the legend (make it smaller) ...
%pos = get(l, 'position');
%pos(3) = 0.75*pos(3);
%set(l, 'position', pos);

xlbl = xlabel('\it Tiempo entre arribos');
ylbl = ylabel('\it Frecuencia');
set(xlbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
set(ylbl, 'Interpreter', 'tex', 'FontName', 'Times', 'FontSize', iStdSize);
axis normal;
grid on;

%legend('Histograma de tiempos entre llegadas al registro');
%xlabel('Tiempo entre arribos');
%ylabel('Frecuencia');

hold off;
%print -deps 'Histograma_E3.eps';


%% Kolmogorov-Smirnov Test:
% matlab function for KS-test ... 
% (to verify the result of manually calculated version) ...
disp('KS-Test:');
%cdf_data = normcdf(times, mu, sigma); % other values with a small diffence between the values below ... 
cdf_data = normcdf(times, mean(times), std(times, 1));
[h, p, ksstat, cv] = kstest(times, [times', cdf_data'], alpha)
% check null hypothesis ...
checkNullHypothesis(ksstat, cv);

% plotting the test statistics of the KS-test ...
figure(3);
clf;
xx = linspace(1.7, 2.3, length(times));
f = cdfplot(times);
hold on;
g = plot(xx, normcdf(xx, mean(times), std(times, 1)), 'k--');
%g = plot(xx, normcdf(xx), 'k--');
set(f, 'Color', 'black', 'LineWidth', 1);
set(g, 'LineWidth', 1);
l = legend([f g], '\it Empirico','\it CDF');
set( l, 'Interpreter', 'tex', 'Location', 'NorthWest', ...
        'FontName', 'Times', 'FontSize', 9 );
% create an annotation ...
x = [0.5 0.4];
y = [0.21 0.26];
txtar = annotation( 'textarrow', x, y, 'HeadWidth', 6, 'interpreter', 'latex', 'string', ...
                    '$\sup_{x}|F_{n}(x) - F(x)|$','FontSize', 11);
hold off;
     
% finally ...
fclose(report);


% scott's rule to calculate the bin width ...
factor = 0.7;
w = 3.49*std(times)*length(times)^(-1/3);
w = w * factor;
nBins = w*length(times)

% without a factor value ...
[w, bins] = scott(times, 1, length(times))



%% result:
% nr. of classes: 7
% df = 7 - 2 = 5 
% 
% j	O(j)	M(j)		X(j)		a		b		E(j)
% 1	7	1.811029	12.677200	1.779700	1.842357	11.934237
% 2	32	1.873686	59.957943	1.842357	1.905014	29.928077
% 3	36	1.936343	69.708343	1.905014	1.967671	48.543521
% 4	52	1.999000	103.948000	1.967671	2.030329	50.945747
% 5	41	2.061657	84.527943	2.030329	2.092986	34.595902
% 6	24	2.124314	50.983543	2.092986	2.155643	15.196929
% 7	6	2.186971	13.121829	2.155643	2.218300	4.315646
% 
% alpha =
% 
%     0.0500
% 
% 
% df =
% 
%      5
% 
% The test does not fit!
% Rejecting the null hypothesis!
% (12.388723 >= 11.070000)
% 
% 
% O =
% 
%      7    32    36    52    41    24     6
% 
% 
% E =
% 
%    11.9342   29.9281   48.5435   50.9457   34.5959   15.1969    4.3156
% 
% 
% h =
% 
%      1
% 
% 
% p =
% 
%     0.0155
% 
% 
% stats = 
% 
%     chi2stat: 12.2687
%           df: 4
%        edges: [1.8110 1.8737 1.9363 1.9990 2.0617 2.1243 2.2496]
%            O: [7 32 36 52 41 30]
%            E: [11.9342 29.9281 48.5435 50.9457 34.5959 19.5126]
% 
% The test does not fit!
% Rejecting the null hypothesis!
% (12.268690 >= 9.490000)
% 
% 
% h2 =
% 
%      0
% 
% 
% p2 =
% 
%     0.1829
% 
% 
% stats2 = 
% 
%     chi2stat: 6.2258
%           df: 4
%        edges: [1.7797 1.8424 1.9050 1.9677 2.0303 2.0930 2.1556 2.2183]
%            O: [8 32 36 52 41 24 7]
%            E: [9.0306 22.5541 43.8987 54.0096 42.0110 20.6553 7.8408]
% 
% The value does not exceed the critical value c (6.225752 < 9.490000).
% Null hypothesis accepted!
% 
% KS-Test:
% 
% h =
% 
%      0
% 
% 
% p =
% 
%     0.7388
% 
% 
% ksstat =
% 
%     0.0475
% 
% 
% cv =
% 
%     0.0952
% 
% The value does not exceed the critical value c (0.047516 < 0.095159).
% Null hypothesis accepted!
% 
% 
% nBins =
% 
%     7.5976
% 
% 
% w =
% 
%    10.8537
% 
% 
% bins =
% 
%    Empty matrix: 1-by-0
